/*
CC=================================================================CC
CC                                                                 CC
CC  Subroutine RSRFFT(X,M):                                        CC
CC      A real-valued, in-place, split-radix FFT program           CC
CC      Decimation-in-time, cos/sin in second loop                 CC
CC      and computed recursively                                   CC
CC      Output in order:                                           CC
CC              [ Re(0),Re(1),....,Re(N/2),Im(N/2-1),...Im(1)]     CC
CC                                                                 CC
CC  Input/output                                                   CC
CC      X    Array of input/output (length >= N)                   CC
CC      M    Transform length is N=2**M                            CC
CC                                                                 CC
CC  Calls:                                                         CC
CC      RSTAGE,RBITREV                                             CC
CC                                                                 CC
CC  Author:                                                        CC
CC      H.V. Sorensen,   University of Pennsylvania,  Oct. 1985    CC
CC                       Arpa address: hvs@ee.upenn.edu            CC
CC  Modified:                                                      CC
CC      F. Bonzanigo,    ETH-Zurich,                  Sep. 1986    CC
CC      H.V. Sorensen,   University of Pennsylvania,  Mar. 1987    CC
CC      H.V. Sorensen,   University of Pennsylvania,  Oct. 1987    CC
CC                                                                 CC
CC  Reference:                                                     CC
CC      Sorensen, Jones, Heideman, Burrus :"Real-valued fast       CC
CC      Fourier transform algorithms", IEEE Tran. ASSP,            CC
CC      Vol. ASSP-35, No. 6, pp. 849-864, June 1987                CC
CC                                                                 CC
CC      This program may be used and distributed freely as         CC
CC      as long as this header is included                         CC
CC                                                                 CC
CC=================================================================CC
*/
#include <math.h>


void rbitrev (double *x, int m);
void rstage (int n, int n2, int n4, double *x1, double *x2, double *x3, double *x4);

void rsrfft (x, m)
double *x;
int m;
{
   int n, is, id, i0, n2, n4, k;
   double t1;
   x--;
   n = 1 << m;
/*C-------Digit reverse counter---------------------------------------C*/
   rbitrev (x, m);
/*C-----Length two butterflies----------------------------------------C*/
   is = 1;
   id = 4;
   do {
      for (i0 = is; i0 <= n; i0 += id) {
         t1 = x[i0];
         x[i0] = t1 + x[i0 + 1];
         x[i0 + 1] = t1 - x[i0 + 1];
      }
      is = 2 * id - 1;
      id = 4 * id;
   } while (is < n);
/*C-------L shaped butterflies----------------------------------------C*/
   n2 = 2;
   for (k = 2; k <= m; k++) {
      n2 = n2 * 2;
      n4 = n2 / 4;
      rstage (n, n2, n4, &x[0], &x[n4], &x[2 * n4], &x[3 * n4]);
   }
   return;
}

/*
C===================================================================C
C  Subroutine RSTAGE - the work-horse of the RFFT                   C
C       Computes a stage of a real-valued split-radix length N      C
C       transform.                                                  C
C  Author                                                           C
C       H.V. Sorensen,   University of Pennsylvania,  Mar. 1987     C
C===================================================================C
*/
void rstage (n, n2, n4, x1, x2, x3, x4)
int n, n2, n4;
double *x1, *x2, *x3, *x4;
{
   int n8, is, id, i1, i2, j, jn;
   double t1, t2, e, t3, t4, t5;
   double ss1, sd1, ss3, sd3;
   double cc1, cd1, cc3, cd3;
   n8 = n2 / 8;
   is = 0;
   id = n2 * 2;

   do {
      for (i1 = is + 1; i1 <= n; i1 += id) {
         t1 = x4[i1] + x3[i1];
         x4[i1] = x4[i1] - x3[i1];
         x3[i1] = x1[i1] - t1;
         x1[i1] = x1[i1] + t1;
      }
      is = 2 * id - n2;
      id = 4 * id;
   } while (is < n);
    /*C*/ if (n4 - 1 <= 0)
      return;

   is = 0;
   id = n2 * 2;
   do {
      for (i2 = is + 1 + n8; i2 <= n; i2 += id) {
         t1 = (x3[i2] + x4[i2]) * .7071067811865475;
         t2 = (x3[i2] - x4[i2]) * .7071067811865475;
         x4[i2] = x2[i2] - t1;
         x3[i2] = -x2[i2] - t1;
         x2[i2] = x1[i2] - t2;
         x1[i2] = x1[i2] + t2;
      }
      is = 2 * id - n2;
      id = 4 * id;
   } while (is < n);
    /*C*/ if (n8 - 1 <= 0)
      return;

   e = 2 * 3.14159265358979323 / n2;
   ss1 = sin (e);
   sd1 = ss1;
   sd3 = 3 * sd1 - 4 * (sd1 * sd1 * sd1);
   ss3 = sd3;
   cc1 = cos (e);
   cd1 = cc1;
   cd3 = 4 * (cd1 * cd1 * cd1) - 3. * cd1;
   cc3 = cd3;

   for (j = 2; j <= n8; j++) {
      is = 0;
      id = 2 * n2;
      jn = n4 - 2 * j + 2;
      do {
         for (i1 = is + j; i1 <= n; i1 += id) {
            i2 = i1 + jn;
            t1 = x3[i1] * cc1 + x3[i2] * ss1;
            t2 = x3[i2] * cc1 - x3[i1] * ss1;
            t3 = x4[i1] * cc3 + x4[i2] * ss3;
            t4 = x4[i2] * cc3 - x4[i1] * ss3;
            t5 = t1 + t3;
            t3 = t1 - t3;
            t1 = t2 + t4;
            t4 = t2 - t4;
            x3[i1] = t1 - x2[i2];
            x4[i2] = t1 + x2[i2];
            x3[i2] = -x2[i1] - t3;
            x4[i1] = x2[i1] - t3;
            x2[i2] = x1[i1] - t5;
            x1[i1] = x1[i1] + t5;
            x2[i1] = x1[i2] + t4;
            x1[i2] = x1[i2] - t4;
         }
         is = 2 * id - n2;
         id = 4 * id;
      } while (is < n);
       /*C*/ t1 = cc1 * cd1 - ss1 * sd1;
      ss1 = cc1 * sd1 + ss1 * cd1;
      cc1 = t1;
      t3 = cc3 * cd3 - ss3 * sd3;
      ss3 = cc3 * sd3 + ss3 * cd3;
      cc3 = t3;
   }
    /*C*/ return;
}


/*
CC=================================================================CC
CC                                                                 CC
CC Subroutine RBITREV(X,M):                                        CC
CC      Bitreverses the array X of length 2**M. It generates a     CC
CC      table ITAB (minimum length is SQRT(2**M) if M is even      CC
CC      or SQRT(2*2**M) if M is odd). ITAB need only be generated  CC
CC      once for a given transform length.                         CC
CC                                                                 CC
CC Author:                                                         CC
CC      H.V. Sorensen,   University of Pennsylvania,  Aug. 1987    CC
CC                       Arpa address: hvs@ee.upenn.edu            CC
CC                                                                 CC
CC      This program may be used and distributed freely as long    CC
CC      as this header is included.                                CC
CC                                                                 CC
CC=================================================================CC
*/
void rbitrev (x, m)
double *x;
int m;
{
   int itab[66000];
   int m2, nbit, imax, lbss, i, j, k, l, j0;
   double t1;
/*C-------Initialization of ITAB array--------------------------------C*/
   m2 = m / 2;
   nbit = 1 << m2;
   if (2 * m2 != m)
      m2 = m2 + 1;
   itab[1] = 0;
   itab[2] = 1;
   imax = 1;

   for (lbss = 2; lbss <= m2; lbss++) {
      imax = 2 * imax;
      for (i = 1; i <= imax; i++) {
         itab[i] = 2 * itab[i];
         itab[i + imax] = 1 + itab[i];
      }
   }

/*C-------The actual bitreversal--------------------------------------C*/
   for (k = 2; k <= nbit; k++) {
      j0 = nbit * itab[k] + 1;
      i = k;
      j = j0;
      for (l = 2; l <= itab[k] + 1; l++) {
         t1 = x[i];
         x[i] = x[j];
         x[j] = t1;
         i = i + nbit;
         j = j0 + itab[l];
      }
   }
   return;
}
